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Abstract. We measure the phase diagram of a 87 Rb Bose gas in a harmonic trap 
in terms of macroscopic parameters obtained from the spatial distribution of atoms. 
Considering the relevant variables as size of the cloud V, number of atoms N and 
temperature T, a novel parameter n = II (AT, V, T) is introduced to characterize 
the overall pressure of the system. We construct the phase diagram (II vs T) 
identifying new features related to Bose-Einstein condensation (BEC) transition in 
a trapped gas. A thermodynamic description of the phase transition based on purely 
macroscopic parameters, provide us with properties that do not need the local density 
approximation. An unexpected consequence of this analysis is the suggestion that BEC 
appears as a continuous third-order phase transition instead of being a second-order 
one. 
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Experiments with cold gases are normally performed in traps which result in 
spatially inhomogeneous samples[U [2j [31 HJ E]- In such cases, the phase diagram for 
the transition from condensate to non-condensate phases becomes difficult to express 
using quantities that apply mostly for homogenous systems, say, diagrams in terms 
of pressure versus temperature. Part of the problem can be solved if one relies on a 
local-density approximation (LDA) where characterization of a non-uniform fluid, at a 
given spatial position, is obtained considering local quantities which in turn determine 
bulk properties. Since the phenomena observed, such as Bose-Einstein condensation 
(BEC), occur at a macroscopic level in a non- uniform trap, the description would be 
more appropriate if the overall characterization of the system could be achieved based on 
few thermodynamic parameters, which in turn should depend on other thermodynamic 
variables such as temperature, number of particles and properties of the confining 
potential. 

Recently, several papers have explored macroscopic quantum phenomena based on 
measurements of thermodynamic properties. Horikoshi and collaborators |6] have used 
a force balance equation to create a state equation, which allow to determine local 
quantities of a Fermi gas at unitarity. Ho and Zhou [7] have proposed an algorithm 
to obtain the phase diagram and thermodynamic quantities of the corresponding bulk 
system, by analyzing the density profile of trapped gases. Equivalent concepts were also 
used by Nascimbene et al.[8] to explore thermodynamics of a Fermi system. 

In this work, we consider a different way to describe the phase diagram for the BEC 
transition, using a single macroscopic parameter to characterize the inhomogeneous 
system. This approach is based on concepts introduced in Refs. [9j [TOj [TTJ 0] , and it 
is applied to experimentally analyze the BEC transition for a harmonic trapped Rb 
gas. However, the employed ideas are quite general and can be applied for any trapped 
system. 

We start by considering the relevant macroscopic parameters to be used. The 
thermodynamic state of a gas trapped in a macroscopic harmonic potential is naturally 
characterized by the trapped number of atoms N, the temperature T and parameters 
that determine the confining potential. In a uniform gas, apart from temperature and 
number of particles, the mechanical variable that specifies the thermodynamic state 
of the system is the volume V that the sample occupies. Since a gas confined in a 
harmonic trap does not have a well defined volume, one can use an alternative variable 
that yields analogous information regarding its size. Considering the harmonic potential 
Vext(f) — ■ r) 2 with Q = (u x ,u y ,u z ), it turns out that in a first approximation 

the size of the system is given by r* 3 [ksT /moj 2 ) 3 ^ 2 . Thus, apart from constants 
and other independent thermodynamic variables, this allows us to identify V = 1/u 3 
as the "volume parameter", where u 3 = u x oj y u> z . This variable plays the role of the 
actual volume in a harmonically trapped fluid. In fact, long ago, de Groot et al.[12j, 
and later Bagnato et al.[13] and others [HI [151 IB], demonstrated that V is the correct 
thermodynamic quantity determined by the external potential, in the sense that the 
thermodynamic limit is obtained when N — » oo, V — > oo but N/V — > constant. The 
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volume parameter, which should not be identified with the real volume of the system, 
does have equivalent properties to the spatial volume of a uniform system. Among the 
most relevant ones is its extensive thermodynamic quality [91 fTU| fTT] . 

Therefore, we take the collection (JV, T, V) as the independent thermodynamic 
variables that specify the state of the system. Any other thermodynamic property 
of the gas can be expressed as a function of those variables. That is, just as in a 
homogeneous fluid the phase diagram may be expressed in terms of the equation of 
state p = p(N, T, V), with p the hydrostatic pressure, here, because the state is given in 
terms of the "volume parameter", we also define a "pressure parameter" IT, such that the 
equation of state becomes now II = T1(N, T, V). Analogously to the volume parameter, 
II is not an actual hydrostatic pressure, although it does retain characteristics of a 
pressure, as we argue now. By considering Helmholtz free energy F as F = F(N, T, V), 
one can introduce II as the variable conjugate to V, namely, 



We note that while II neither has units of pressure nor V units of volume, the product 
IIV still has units of energy. Clearly, II is an intensive thermodynamic variable and, 
therefore, its dependence on (N,T,V) must be of the form, II = U(N/V,T). Thus, in 
keeping with the analogy, we shall call N/V the density parameter. From the above 
definition of II, without resorting to LDA, a very useful and general expression can be 
derived, this is® EH Q3], 

2 f 1 

II = — / n(r) ~m(uj ■ r) 2 d 3 r, (2) 



3V J w 2 

where n{f) is the inhomogeneous density profile of the gas. Since both n(r) and the 
confining potential are measurable quantities, II can be readily determined. With the 
equation of state IT = U(N/V,T), the phase diagram for the transition can be fully 
described based on measured macroscopic quantities of the global system, without the 
necessity of considering the corresponding local uniform system as done with LDA. 

An additional insight into the physics of II can be gained by noticing the additional 
identity PUJ, 

nv = y p(f)d\ (3) 

showing that, actually, the pressure parameter II is proportional to the integral of the 
local hydrostatic pressure p(r) over all space. To obtain this result we just need to 
remember that mechanical equilibrium in the fluid is determined by the balance of local 
forces, namely, Vp(r) + n{f)W ext {f) = 0. By using equation (j2J) the integration of the 
virial of these forces, r ■ f, leads to (J3J). A similar mechanical equilibrium argument has 
also been used by Horikoshi et al. [6]. 

To measure the quantities involved in Eq.(j2]) across BEC transition we have used 
an experimental setup composed of a double MOT and a QUIC-type trap. The system 
has been described in detail elsewhere |1 7] . In brief, we collect 10 9 87 Rb atoms in a 
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MOT, producing a sample around 100/iK. After a population transfer to the |2, 2) 
hyperfine state, the sample is used to load a QUIC magnetic trap. While trapped, a RF- 
induced evaporation is applied to obtain quantum degeneracy. The trapping potential 
is harmonic with frequencies u x = Uq and u y = u z = 9uq, where u>o = 2tt x 23 Hz. The 
characterization of the sample is done by absorption imaging it on a CCD camera after 
15 ms time of flight. We typically produce BEC samples containing 1-7 xlO 5 atoms. 

Each acquired image is fitted with a bimodal distribution composed of a gaussian 
and an inverted parabola (Thomas- Fermi distribution) which represent, respectively, 
the thermal cloud and the condensate part. By using the procedure described by Castin 
and DumpS], one can obtain the "in situ" density distribution of the cloud, n(f), which 
is cylindrically symmetric. From those fittings, the temperature T is obtained (gaussian 
tail), as well as the condensate fraction Nq and the total number of atoms N through a 
spatial integration. Finally, by using equation (j2J) we obtain the pressure parameter IT. 

For the experiment here reported we have kept the trap frequencies fixed, that 
is V = constant. Temperature and number of atoms were varied by controlling the 
evaporative cooling final frequency and the time of evaporation [17] . We have performed 
about 300 runs of the experiments and arranged the results with similar number of 
particles for different temperatures. 

Since we are expressing II = U(N/V, T), the experiment corresponds to the 
investigation of II vs T, while the density parameter N/V remains constant. The results 
of this experiment is shown in Fig. [TJ For a fixed density parameter we observe two 
different regions of behavior as T decreases. First, for high temperatures, II has an 
almost linear dependence both on T and N/V (which is the prediction of an ideal 
Bose gas above BEC), until an abrupt change takes place, and the decrease of II 
with T becomes faster than linear. This change of behavior characterizes the critical 
temperture T c = T C (N/V) where BEC takes place. This change was also detected in the 
corresponding density profiles, that is, those below T c showed the characteristic peak of 
the condensate. The critical point (Il c , T c ) for each density parameter was determined 
by the crossing point of the extrapolations of both sides of the transition. Following the 
definition of II from equation ([TJ, and the fact that for T > T c the density clouds are very 
well fitted by a gaussian profile, the mentioned linear dependence is expected. That is, 
if one considers the Bose distribution of particles in the non-interacting approximation, 



with g n (z) the Bose function and z = exp(/i/fceT) the fugacity, the deviation from 
linearity in T and N/V is rather small since 1 > g^z) / g 3 (z) > 0.9. Even if Hartree- 
Fock corrections are included [PJ, the dependence of II on T remains very close to linear 
for T > T c . Below T c , the behavior of II is intrinsically related to the superfluid nature of 
the gas, depending very strongly now on the magnitude of the interparticle interaction. 
That is, if the interparticle interactions played no role, all the curves in Fig. [TJ would 
follow the same ideal transition line (similar to the dash line in Fig. [21 see below). 

From the critical points determined in Fig. [Tj one can extract the phase diagram IT 



n 



Nk B T gi (z) 
V g 3 (z) 



(ideal gas), 



(4) 




Figure 1. Pressure parameter II as a function of temperature T for different 
constant number of particles N, crossing their corresponding critical temperature 
T C (N/V). Solid lines are best fittings through experimental points. Above the critical 
temperature the dependence is linear. 



vs T shown in Fig. El where the critical line Il c vs T c separates the thermodynamic states 
into two domains, one where the fluid is thermal overall, and the other where there is a 
superposition of a Bose condensate fluid with a thermal component. The line U c vs T c 
represents what is expected to be a continuous phase transition[16j between a normal 
gas and a superfluid in an interacting Bose fluid confined by a harmonic trap. Below 
we address the issue of the order of the transition as predicted by Fig. [U If the gas 
were ideal one should obtain Il c ~ T c 4 at the transition line, however, a logarithm plot 
reveals that a simple law I1T 7 = constant is not obeyed for 7 constant. The exponent 
becomes larger as the temperature decreases. This feature, and those discussed below, 
are certainly due to the influence of interactions below T c and to the presence of the 
confining potential. 

We now analyze data for each isodensity curve (density parameter constant) shown 
in Fig. (H for temperatures below T c . These points are presented in Fig. [31 Note 
that higher isodensity curves produce higher pressure parameter values at a given 
temperature. This appears as expected from an overall analysis of these variables in 
harmonic [9] and linear quadrupolarjll] traps using the Hartree-Fock approximation. 
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Figure 2. Phase diagram in terms of pressure parameter II and temperature T. The 
BEC transition line separates the two phases: purely thermal cloud and thermal cloud 
plus condensate. 



The extrapolation of the experimental points in Fig. [31 represented by the solid lines, 
allows the determination of il = U(N/V,T = 0), the zero-temperature pressure 
parameter, which depends on the pure condensate density without the presence of the 
thermal component. At T = 0, the Thomas- Fermi (TF) approximation works on its 
best conditions. Therefore, by using TF and LDA|14l [15] and assuming a cloud of iV 
interacting atoms with scattering length a, equation ()2]) at T = yields, 

n = i (i5h 2 aV^) 2/5 (Nu 2 r co z ) 7/5 . (5) 

It is found then, that at T = the pressure parameter scales with the number of particles 
(keeping V fixed) as n ~ N s . A log-log plot of the values found from Fig. [3] yields 
the value 5 = 1.5lg;2, we ^ within the expected value 1.4 predicted by TF. The use of 
expression ([5]) together with the data presented in Fig. El allow the determination of 
the scattering length, fully based on global thermodynamic quantities. However, this is 
not the scope of this paper. 

As with standard pressure and volume in uniform fluid systems, the variables II and 
V provide us with a macroscopic thermodynamic description of the fluid without the 
necessity of relying on LDA. Because of the physical connection between the parameters 
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Figure 3. Detail of pressure parameter dependence on temperature, for different 
fixed number of particles, If = H(N/V,T), in the BEC region. The extrapolation 
of the experimental points allows for the determination of the pressure parameter at 
T = 0. 

(IT, V) with (p, V), the phase diagram reveals interesting macroscopic features not yet 
considered along the phase transition for the formation of a superffuid atomic phase in 
a trap. 

As a first illustration, we note the striking similarity of the phase diagram of Fig. [2] 
with the normal gas to superfluid phase transition in the bottom part of the phase 
diagram p — Toi 4 He, see Refs. [T9l 120] for instance. We point out that this transition is 
different from the A-transition which occurs between a normal liquid and a superfluid. In 
the gas-superfluid transition here seen the coexistence pressure parameter vanishes as the 
temperature goes to zero, just as in the analogous gas-superfluid transition diagram p—T 
of 4 He. On the other hand, we also insist in the profound differences between considering 
a local picture and a global or thermodynamic one of a confined system. If studies were 
made of local variables only, the phase diagram obtained would refer to the bulk system 
and the transition line would not yield the information of the occurrence of BEC in the 
trapped gas. That is, knowledge of the bulk phase diagram solely does not allow a direct 
calculation of the phase diagram of the inhomogeneous trapped gas. One needs first 
full knowledge of the number particle density n = N/V and hydrostatic pressure p of 
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the homogeneous gas as a function of both chemical potential /x and temperature T, i.e. 
n(/i, T) and p(fi,T). Then, using LDA with the harmonic potential, a reconstruction 
of the phase diagram n vs T such as in Fig. [1] could be made. However, BEC occurs 
in the trapped inhomogeneous gas and this information is succinctly contained in the 
phase diagram of Fig. CD and Fig. EJ 

The elucidation of the order of a phase transition is a very delicate issue that can 
only be resolved by either solving the problem exactly or by resorting to experimental 
evidence. The liquid to superffuid transition in 4 He is a continuous second order 
one [21] and the BEC transition in bulk systems is widely believed to be also of second 
order [22} [16]. The features that determine this behavior is the divergence of the heat 
capacity and of the isothermal compressibility [23]. The former is well exemplified by 
the A-transition in 4 He. Since the isothermal compressibility is related to the density 
or particle fluctuations, its divergence signals the appearance of critical behavior in 
the system [23]. Here, by having knowledge of the equation of state n = U(N/V,T) 
through the plots of Fig. CD, we can attempt to calculate the corresponding isothermal 
compressibility parameter K T , which turns out to be related to the density fluctuations 
as well, that is, 



T N\dN/VJr 

V 1 AiV 2 



Nk B T N 



(6) 



-2 



where the particle fluctuations are AN 2 = N 2 — iV . The above expressions follow 
from the definition of n, equation ([T]). From the data of Fig. [1] we were able to 
obtain the isothermal compressibility for one isodensity curve (N ~ 2.8 x 10 5 ) as 
a function of temperature. The result is shown in Fig. HI where we also show the 
approximated location of the corresponding BEC temperature T c . Two features must 
be highlighted. First, we see no evidence that the isothermal compressibility parameter 
becomes very large at criticality (diverging in the thermodynamic limit) neither from 
above nor from below T c , as expected in a second-order phase transition [23]. But 
secondly and very interestingly, one sees a change of sign in the curvature of Kt as 
a function of T, the zero curvature point nearly coinciding with the location of T c . 
A plot of the derivative {dKrp/dT)^/y will then show a pronounced (negative) peak 
at T c . Although it is widely believed that phase transitions can only be of first or 
second order, this knowledge is solely based on uniform systems. Here, we advance the 
suggestion that BEC transition in a harmonically trapped gas may be of third order. We 
make this possible identification following Ehrenfest classification. That is, we identify 
that derivatives of the compressibility can be expressed as third-order derivatives of a 
free energy, 

(dK T \ 1 / d 3 G /dU\ d 3 G 

\&T) N/V ~ "V [WdlP + \fff) m dn* 

where G = G(N,H,T) is Gibbs free energy. If the peak of (dK T /dT) N / v at T c 



(7) 
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becomes more and more pronounced, as the system becomes larger (i.e. approaching 
the thermodynamic limit), it would indicate the divergence of a third derivative of the 
free energy. This would also indicate that while the density fluctuations do not diverge, 
their derivatives do so. This, to the best of our knowledge, is a property never reported 
before. We believe that this behavior is due to the fact that the system is confined by a 
harmonic potential and that it is not uniform. This behavior can be seen to be already 
present as a precursor in the ideal Bose gas in a harmonic trap. That is, one can show 
that for T > T c , it is obeyed, 

/ r)K ideal \ /V\ 2 k n T 

Approaching BEC from above for N/V fixed, i.e. z — > 1 + , one finds that both, #3(1) = 
C(3) and g 2 {l) = vr 2 /6, are finite, while #i(l) ->■ 00. Thus {dK l T deal /dT) N/v -> -00 as 
the experiments suggest. The full elucidation of this aspect requires a much thorough 
analysis of the equation of state and will be reported elsewhere. 

The capability to describe an inhomogeneous quantum gas via macroscopic 
parameters creates means to access other thermodynamic quantities proper of the 
trapped gas that cannot be accessed using local variables within LDA. For instance, 
the heat capacity at constant V, Cy, is a quantity that pertains to the full macroscopic 
system and cannot be directly obtained from direct integration of local quantities [TT]. It 
can be shown that the heat capacity may be obtained by full knowledge of the equation 
of state n = U(N/V,T) and by independent measurements of the adiabatic change 
of the temperature as the volume parameter is changed Further, as done in this 
paper, one can also calculate the corresponding isothermal compressibility parameter 
K T = (N /V)(dH/ 'd(N /V))t • The elucidation of these quantities, K T and Cy, through 
BEC, together with the phase diagram here shown, would reveal the critical details of 
interacting BEC in a trap, a question still unanswered to this day. 
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